A previous presentation developed an approximation to the orbits of a power potential analytic in the exponent of the potential using elliptic functions. While quite accurate for positive exponents, the approximation is visibly deficient for negative exponents. An alternate approach in this case is an approximation form introduced by Struck and extended by Lynden-Bell. One major difference in taking this approach is that temporal dependence is not available, since the orbits are parametrized by angle rather than time.

First some orientation: the Hamiltonian for a body moving under the influence of a spherically symmetric n-dimensional power potential with an arbitrary real exponent is

$H=\frac{p^2}{2m}+a{r}^k=E$

where the variables have their usual definitions as vector squares:

$r^2=(\mathbf{r}·\mathbf{r})=\sum _{i=1}^n{r}_i^2{p}^2=(\mathbf{p}·\mathbf{p})=\sum _{i=1}^n{p}_i^2$

Since the motion takes place in an invariant plane of two dimensions, the overall dimensionality of the space will not enter the development. Replacing momenta conjugate to angular variables with the square of angular momentum defines an energy equation for the radial variable:

$\frac{1}{2}m{\stackrel{·}{r}}^2+\frac{L^2}{2m{r}^2}+a{r}^k=E$

This equation can be rearranged, integrated and in principle inverted to produce the radial variable as a function of time, but what is sought here is a parametrization by angle. Angular momentum conservation

$m{r}^2\stackrel{·}{\phi }=L$

provides a differential equation linking a polar angular variable and the radial variable. This can be used to rewrite the energy equation as

$\frac{dr}{d\phi }=\frac{r}{L}\sqrt{2mE{r}^2-2ma{r}^{k+2}-L^2}$

This equation is not directly suitable for integration due to the indeterminancy of the sign of the square root. Instead take another derivative with respect to angle and use this equation to replace all first derivatives:

$\frac{d^{2}r}{d{\phi }^2}=\frac{4mE}{L^2}{r}^3-\frac{ma(k+4)}{L^2}{r}^{k+3}-r$

A numerical integration of this equation for k > 0 with typical initial conditions and unit mass looks like this parametrically:

For k < 0 the coupling constant must also be negative for bound orbits. With an angular momentum dependent on exponent, typical parametric orbits look like this:

The orbits in these interactive graphics start at a point farthest from the center of attraction for numerical convenience. With both graphics one can see a variety of closed orbits for nontrivial exponents, which is discussed in another presentation.

Following Lynden-Bell with adjustments for consistency with previous presentations, the analytic approximation to be considered is of the form

$(\frac{l}{r}{)}^{k+2}=1-ecosn\phi$

The generalized semilatus rectum l and the generalized eccentricity e are given in terms of the points of furthest and nearest approach to the center of attraction, designated here as r₀ and r₁. The expressions given by Lynden-Bell are

$l=[\frac{1}{2r_0^{k+2}}+\frac{1}{2r_1^{k+2}}{]}^{-\frac{1}{k+2}}e=\frac{r_0^{k+2}-r_1^{k+2}}{r_0^{k+2}+r_1^{k+2}}$

The approximation is identically satisfied at r₀ with these expressions and $\phi =0$ .

The points r₀ and r₁ are extrema at which the radial motion has a vanishing temporal derivative. Given either extremum and an angular momentum L, the energy equation can be used to determine the value of the energy and then the location of the other extremum:

$\frac{L^2}{2m{r}_0^2}+a{r}_0^k=E=\frac{L^2}{2m{r}_1^2}+a{r}_1^k$

While it is possible to develop a power series expression for the second extremum, in practice numerical root finding is more convenient and accurate.

All that remains is to determine the remaining parameter n in the approximation. Struck and Lynden-Bell present various ways to do this with approximate solutions to the underlying differential equation for the radial variable. This parameter will be chosen here to make the orbit approximation match the actual motion at the extrema.

To do so one first needs the angular change in moving between the two extrema. The angular position of the orbiting body is a simple result of angular momentum conservation:

$m{r}^2\stackrel{·}{\phi }=L\to \phi =\int dt\frac{L}{m{r}^2}$

The integral can be written entirely in terms of the radial variable using the energy equation. The angular change required is

$\Delta \phi =\underset{r_1}{\overset{r_0}{\int }}\frac{Ldr}{r\sqrt{2mE{r}^2-2ma{r}^{k+2}-L^2}}$

where the order of integration is chosen to keep the result positive. Putting this angular value into the analytic approximation at r₁ gives

$n=\frac{1}{\Delta \phi }{cos}^{-1}\frac{1}{e}[1-(\frac{l}{r_1}{)}^{k+2}]$

and all of the pieces are in place.

For k < 0 the approximation compared with the numerical integration is

The input parameters are limited in range to avoid numerical errors. There is enough range in the exponent of the potential to see clearly the transition from the exact ellipse of the Kepler potential to rotating orbits. Variation of the angular momentum shows, as pointed out by Lynden-Bell, that the approximation is most accurate for lower eccentricity, or higher angular momentum. For higher eccentricity, or lower angular momentum, the approximation is noticeably different from the numerical integration.

Naïvely attempting to use the same approximation for k > 0 results in behavior that clearly does not match the numerical integration except for high angular momentum:

This is because the approximation must assume a different form under a transformation group that takes power potentials with negative exponents over into positive exponents. Consider a new set of variables defined by

$\frac{d\overline{r}}{d\overline{t}}=f\left(r\right)\frac{dr}{dt}$

The energy equation without an explicit potential is in these variables

$\frac{m}{2}(\frac{d\overline{r}}{d\overline{t}}{)}^2=\frac{m}{2}{f}^2(\frac{dr}{dt}{)}^2=f^2[E-\Phi \left(r\right)-\frac{L^2}{2m{r}^2}]\equiv \overline{E}-\overline{\Phi }\left(\overline{r}\right)-\frac{L^2}{2m{\overline{r}}^2}$

Equating energy terms with potential terms, as well as equating angular momentum terms, gives the transformation

$f^2=-\frac{\overline{\Phi }\left(\overline{r}\right)}{E}=-\frac{\overline{E}}{\Phi \left(r\right)}=\frac{r^2}{{\overline{r}}^2}$

For the power potentials in this presentation, the transformation of the radial variable becomes

$\overline{r}=\sqrt{-\frac{\alpha }{\overline{E}}}{r}^{1+k/2}$

This expression remains real because energy and coupling constant have the same sign for bound orbits in the same variables, but opposite signs in combinations of barred and unbarred variables. The transformation for the potential itself is

$\begin{array}{l}\overline{\Phi }\left(\overline{r}\right)=-E\frac{r^2}{{\overline{r}}^2}=-\frac{E}{{\overline{r}}^2}(-\frac{\overline{E}}{\alpha }{\overline{r}}^2{)}^{\frac{2}{k+2}}\\ =-E(-\frac{\overline{E}}{\alpha }{)}^{\frac{2}{k+2}}{\overline{r}}^{-2k/(k+2)}\equiv \overline{\alpha }{\overline{r}}^{\overline{k}}\end{array}$

The exponent in this transformation as well as that of the radial variable are identical to those between pairs of power potentials in the Schrödinger equation. Applying the transformation a second time to the radial variable, starting from barred quantities, gives

$\overline{\overline{r}}=\sqrt{-\frac{\overline{\alpha }}{E}}{\overline{r}}^{1+\overline{k}/2}=(-\frac{\overline{E}}{\alpha }{)}^{\frac{1}{k+2}}(\sqrt{-\frac{\alpha }{\overline{E}}}{r}^{1+k/2}{)}^{\frac{2}{k+2}}=r$

Since the transformation returns identically to itself, there does not appear to be a simple way to determine the transformed energy. Fortunately this is not needed when transforming the approximation. The factor on the radial variable containing the transformed energy falls explicitly out of the generalized eccentricity. While it appears in the generalized semilatus rectum, it is of the same power as the transformed radial variable and so will drop out of their ratio. Thus one can simply take $r={\overline{r}}^{2/(k+2)}$ and define the transformed parameters

$\overline{l}=[\frac{1}{2{\overline{r}}_0^2}+\frac{1}{2{\overline{r}}_1^2}{]}^{-1}\overline{e}=\frac{{\overline{r}}_0^2-{\overline{r}}_1^2}{{\overline{r}}_0^2+{\overline{r}}_1^2}$

and since the exponent dependent on k in the original approximation cancels exponents in the the transformed quantities, the approximation now takes the form

$\frac{\overline{l}}{{\overline{r}}^2}=1-\overline{e}cosn\phi$

with the constant exponent on the left-hand side noted by Lynden-Bell. It now remains to consider transformation of the angular variable using the statement of angular momentum conservation:

$d\overline{\phi }=\frac{L}{m{\overline{r}}^2}d\overline{t}=\frac{L}{m{\overline{r}}^2}\times \frac{\overline{r}d\overline{r}}{rdr}dt=\frac{dln\overline{r}}{dlnr}d\phi$

For power potentials this general transformation is simply

$d\overline{\phi }=(1+\frac{k}{2})d\phi =(1-\frac{\overline{k}}{\overline{k}+2})d\phi =\frac{2}{\overline{k}+2}d\phi$

Since the multiplicative factor is constant for a given potential, it can be absorbed into the free parameter n for the purposes of numerical presentation. Evaluation of the angular change between extrema and this final parameter proceed as before, giving a more accurate approximation for positive exponents:

The differential equation defining the exact orbits is numerically unstable around the inner extremum r₁ for low angular momentum, so to avoid blatantly incorrect behavior a lower limit is placed on the angular momentum input. This still allows enough range of initial conditions to get a feel for how well this approximation works for positive exponents.

The approximation is again most accurate for lower eccentricity, or higher angular momentum, while noticeably different from the numerical integration for higher eccentricity, or lower angular momentum. Cursory visual inspection appears to indicate that positive exponents are better approximated by the form using elliptic functions, but this would need to be made numerically explicit.

Having visually presented the difference in approximations for positive exponents before and after the application of the transformation group suggests that a similar transformation would benefit the form using elliptic functions. That will be left for a future presentation.

Uploaded 2019.04.28 analyticphysics.com